function [] = computeCorrelation_TheoryFix( dirName ,fileNum ,ncM)
%COMPUTECORRELATIONOFR_THEORY Summary of this function goes here
%   Detailed explanation goes here
%---------------------------Preparation-------------------------------%
% this part repeats the "Preparation" of plotEntropyFix

fileArray = dir(dirName);
nameCat = [dirName,'\',fileArray(fileNum+3).name];
% nameCat = getName( name ,minute, t ); % for a fixed frame
[position,velocity,num]=readText( nameCat );

S = sum(velocity) / num; % now S is the average velocity vector 
SL = norm(S); % SL equals the model of S
unitVec = S/SL;
sL = velocity * unitVec'; % column vector
pai = velocity - sL * unitVec ;  % get the longitudinal velocity "sL*unitVec" and the perpendicular velocity "pai" 

distance = zeros( num,num );
expCorr = zeros(num,num);
for i= 1 : num
    for j = 1: num
        distance(i,j) = norm(position(i,:)-position(j,:));
        expCorr(i,j) = sum(velocity(i,:).*velocity(j,:));
    end
end % acquire distance and correlation between i and j

[logicI,logicB] = distinguishBoundary(position,distance);
numI = sum(logicI);
distance_II = distance(logicI,logicI);


%-----------------------Compute CorrP--------------------%
deltaR = 0.5;
maxR  =  floor(max(max(distance)))+1;
maxN = floor(maxR/deltaR);

[logic, ~, intCorrelationI] = computeIntCorrelation( num ,ncM ,distance ,expCorr ,logicI);
N_OO = 0.5 * (logic + logic'); % when i = j ,logic(i,j) = False
N_IB = N_OO(logicI,logicB);
N_II = N_OO(logicI,logicI);
AdJ_II = diag(sum(N_II,2) + N_IB * sL(logicB)) - N_II;
%h_I = N_IB * velocity(logicB,:); % column vector
hL_I = N_IB * sL(logicB);
hP_I = N_IB * pai(logicB,:);
J =(numI-1)*( 0.5 *sum(sum( AdJ_II^-1 .*( hP_I * hP_I'))) - 0.5 * (norm(sum(pai(logicB,:))+sum( AdJ_II^-1 * hP_I)))^2 / (sum(sum(AdJ_II^-1))) +...
    sum(hL_I) + 0.5 * sum(sum(N_OO(logicB,logicB).*(velocity(logicB,:) * (velocity(logicB,:))'))) + 0.5 * ncM * (numI - num*intCorrelationI) )^(-1);
disp(J);
%J = 20;
A_II = AdJ_II * J;
hP_I = hP_I * J;

CorrP_TheoR =zeros(1,maxN);
CorrL_TheoR =zeros(1,maxN);
count2 = zeros(1,maxN);
% AA = repmat(sum(A_II^-1,2),1,3).*hP_I;
AA = A_II^-1 * hP_I; % a controversial expression
pai_Theo = AA - sum(A_II^-1,2) * ( sum(pai(logicB,:)) + sum(A_II^-1 * hP_I) ) / sum(sum(A_II^-1)) ;

CorrP_Theo = pai_Theo * pai_Theo' + 2 * ( A_II^-1 - sum(A_II^-1,2) * sum(A_II^-1) / sum(sum(A_II^-1)) ) ;
% CorrP_Theo = pai_Theo * pai_Theo' + 2 * ( A_II^-1 - sum(A_II^-1,2) * sum(A_II^-1) / sum(sum(A_II^-1)) ) ;
CorrL_Theo = (1-0.5*diag(CorrP_Theo)) * (1-0.5*diag(CorrP_Theo))';

for i = 1 : numI
    for j = 1 : numI
        temp = floor((distance_II(i,j)/deltaR))+1;
        if temp < maxN
        CorrP_TheoR(temp) = CorrP_TheoR(temp) + CorrP_Theo(i,j);
        CorrL_TheoR(temp) = CorrL_TheoR(temp) + CorrL_Theo(i,j);
        count2(floor((distance_II(i,j)/deltaR))+1) = count2(floor((distance_II(i,j)/deltaR))+1)+1;
        end
    end
end % transfer (i,j) 's relation to R 's relation

meanCorrP_TheoR = CorrP_TheoR./count2;
meanCorrL_TheoR = CorrL_TheoR./count2;
meanCorrP_TheoR(1) =[]; 
meanCorrL_TheoR(1) =[]; 
meanCorr_TheoR = meanCorrP_TheoR + meanCorrL_TheoR;
%meanCorrP_TheoR(end) = []; % delete the terminal data of min and max
%-------------------------------Plot--------------------------------%

figure(1)
plot(deltaR*(2:maxN),meanCorrL_TheoR,'kx')
figure(2)
plot(deltaR*(2:maxN),meanCorrP_TheoR,'kx')
figure(3)
plot(deltaR*(2:maxN),meanCorr_TheoR,'kx')
end

